Bare Necessities

# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
library(nlme)
library(ggimage)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "white", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom"))



# labels for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Flag effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * area * season) ) %>% 
    mutate(
      area = factor(group, levels = area_levels_long_all),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

Everything in this markdown uses the size spectra exponent values where I shifted the wmin parameter to 16g. The wmin values are fixed at 16g across all groups.

# # Data used for Wigley estimates
# wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>%   
#   left_join(area_df)
# 
#### 1g Cutoff Data  ####
# # Load the three regional datasets
# wigley_bmspectra <- read_csv(
#   here::here("Data/processed/wigley_species_bodymass_spectra.csv")) %>% 
#   mutate(survey_area = factor(survey_area, levels = area_levels))
# 
# 
# # Load Shelfwide
# wigley_bmspectra_shelf <- read_csv(
#   here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))
# 
# 
# 
# # Join them together
# wigley_bmspectra_all <- bind_rows(wigley_bmspectra, wigley_bmspectra_shelf) %>% 
#   mutate(metric = "bodymass_spectra",
#          community = "Wigley Subset",
#          survey_area = factor(survey_area, levels = area_levels_all)) %>% 
#   left_join(area_df) %>% 
#   mutate(area = factor(area, area_levels_long_all))


# Modeling datasets

# 
# # # Model Datasets for Wigley Species Subset
# wigley_bmspectra_model_df <- read_csv(here::here("Data/model_ready/wigley_community_bmspectra_mod.csv"))
# 
# # load shelfwide stuff too
# wigley_bmspectra_shelf <- read_csv(
#   here::here("Data/processed/shelfwide_wigley_species_bodymass_spectra.csv"))
# 
# # Bottom Temperatures
# bot_temps <- read_csv(here::here("Data/processed", "trawl_region_seasonal_bottom_temps.csv"))
# 
# 
# # Bolt that on
# wigley_bmspectra_model_df <- bind_rows(
#   wigley_bmspectra_model_df,
#   left_join(
#     wigley_bmspectra_shelf, 
#     bot_temps, 
#     by  = join_by("est_year" == "year", survey_area, season))) %>% 
#   left_join(area_df) %>% 
#   mutate(area = if_else(survey_area == "Northeast Shelf", "Northeast Shelf", area),
#          area = factor(area, levels = area_levels_long_all),
#   season = fct_rev(season))

# Data used for Wigley estimates
wigley_in <- read_csv(here::here("Data/processed/wigley_species_trawl_data.csv")) %>%   
  left_join(area_df) 

# Data used for 16g wmin
wigley_in_16 <- wigley_in %>% 
  filter(wmin_g >= 16)

# Load the 16g minimum size cutoff
wigley_bmspectra_model_df <- read_csv(
  here::here("Data/model_ready/wigley_community_16g_bmspectra_mod.csv")) %>% 
  mutate(
    metric = "bodymass_spectra",
    community = "Wigley Subset",
    survey_area = factor(survey_area, levels = area_levels_all),
    area = factor(area, area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))
  )

Part 1: Summary of Trends

For the first section the focus is on understanding changes in the size structure. This will be done using three metrics:

1. bodymass spectra
   
2. Size quantiles (10th, median, 90th)
   
3. large and small fish indices

This section should address:

Has the size distribution changed over time, where has it changed, and some nuance around what aspects about the distribution have shifted.

For the figures in part 1: non-zero trends over time have been highlighted in the appropriate plot panels with linear trends overlaid.

A. Size Spectra Trends

Size spectra are steepening in 3 seasonal communities. Gulf of Maine (Spring), Georges Bank (Spring), and Southern New England (Spring).

The size distribution in the Mid-Atlantic Bight (Spring) is becoming more shallow.

# Model linear trends in time
spectra_trend_mod <- gls(
  b ~ scale(est_year) * area * season,
  correlation = corAR1(form = ~ est_year | area/season),
  data = wigley_bmspectra_model_df )

#check_model(spectra_trend_mod)
# plot(check_collinearity(spectra_trend_mod))


# Get the predictions and flag whether they are significant
bmspectra_trend_predictions <- get_preds_and_trendsignif(spectra_trend_mod) %>% 
  mutate(model_id = "bodymass_spectra-Wigley Subset") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric),
         area = factor(area, levels = area_levels_long_all))


# Make the plot
bmspectra_trend_p <- ggplot() +
   geom_ribbon(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    data = wigley_bmspectra_model_df ,
    aes(est_year, b, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = wigley_bmspectra_model_df ,
    aes(est_year, b, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries" ) +
  geom_line(
    data = filter(bmspectra_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(area ~ ., scales = "free") +
  labs(
    x = "Year",
    y = "Size Distribution Exponent",
    linetype = "Trend",
    color = "Season",
    fill = "Season")
 
# Show plot
bmspectra_trend_p

B. Median Weight Trends

Median body weight is declining in 3 seasonal communities, Gulf of Maine (Spring and Fall) and Georges Bank (Spring).

Median weight in the Mid-Atlantic Bight has been increasing in the Spring.

# Load the median weight/length data
wigley_size_df <- read_csv(here::here("Data/processed/wigley_species_size_summary.csv"))

# shelf-wide summary
wigley_sizes_shelf <- read_csv(here::here("Data/processed/shelfwide_wigley_species_size_summary.csv"))


# Join them
size_results <- bind_rows(wigley_size_df, wigley_sizes_shelf) %>% 
  mutate(
    community = "Wigley Subset",
    survey_area = factor(survey_area, levels = area_levels_all)) %>% 
  left_join(area_df)



# Get trends:
wt_mod_wig <- gls(
  med_wt_kg ~ scale(est_year) * area * season,
  correlation = corAR1(form = ~ est_year | area/season),
  data = drop_na(size_results, med_wt_kg))



# Get the predictions and flag whether they are significant
size_trend_predictions <- get_preds_and_trendsignif(wt_mod_wig) %>% 
  mutate(model_id = "med_wt_kg-Wigley Subset") %>% 
  separate(model_id, into = c("metric", "community"), sep = "-") %>% 
  mutate(metric = fct_rev(metric))




# Median length - all species
size_long <- size_results %>% 
  pivot_longer(cols = c(med_len_cm, med_wt_kg), 
               names_to = "metric", 
               values_to = "val")  %>% 
  mutate(
    survey_area = fct_relevel(survey_area, area_levels_all),
    area = fct_relevel(area, area_levels_long_all),
    season = fct_rev(season))


# Median weight - wigley species

# weight plots
wt_plot <- size_long %>% 
  filter(metric == "med_wt_kg") %>% 
  drop_na(val) %>% 
  ggplot() +
  geom_ribbon(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, val, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, val, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_line(
    data = filter(size_trend_predictions, metric == "med_wt_kg", non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_fill_gmri() +
  scale_color_gmri() +
  facet_grid(area ~ .,
             scales = "free")  +
  labs(
    x = "Year",
    y = "Median Weight (kg)",
    color = "Season",
    linetype = "Trend",
    fill = "Season")

wt_plot

Notes:

There are a number of years when median weight surges nearly 0.5kg in the Mid-Atlantic Bight. For these situations I think we likely could trace these down to exceptionally high catches of larger species like elasmobranchs.

We are not seeing a see a shelf-wide decline in size so I might need to check what Friedland did or prepare for a confrontation.

C. Bottom Temperature Trends

Bottom temperature data is from the DuPontavice and Saba combined bottom temperature product. This data product has been clipped to the study areas and seasonally averaged to match the survey sampling seasons.

Plotted below are the long-term trends in each area for each survey season. Significant seasonal trends are shown with linear trends and error displayed.

Seasonal bottom temperatures have been increasing in all 8/10 seasonal communities.

Spring temperature in SNE + MAB have not increased over this period.

# Model linear trends in time
btemp_mods <- lm(
  formula = bot_temp ~ scale(est_year) * area * season,
  data = wigley_bmspectra_model_df)



# change the stupid labels around because
btemp_preds <- as.data.frame(
    # Model predictions
    ggpredict(
        btemp_mods, 
        terms = ~ est_year * area * season) ) %>% 
    mutate(
        area = factor(group, levels = area_levels_long_all),
        season = factor(facet, levels = c("Spring", "Fall")))
btemp_signif <- emtrends(
      object = btemp_mods,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
btemp_trend_predictions <- left_join(btemp_preds, btemp_signif)  %>% 
  mutate(area = group,
         area = factor(area, levels = area_levels_long_all),
         season = factor(season, levels = c("Spring", "Fall")))




# Plot temperatures
temp1 <- wigley_bmspectra_model_df %>% 
  ggplot() +
  geom_ribbon(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, ymin = conf.low, ymax = conf.high, fill = season),
    linewidth = 1, alpha = 0.35) +
  geom_point(
    aes(est_year, bot_temp, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, bot_temp, color = season, linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_line(
    data = filter(btemp_trend_predictions, non_zero == TRUE),
    aes(x, predicted, color = season, linetype = "Significant Trend"),
    linewidth = 1, alpha = 1) +
  scale_color_gmri() + 
  scale_fill_gmri() + 
  #scale_linetype_manual(values = c(1,1)) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = label_number(suffix = "\u00b0C")) +
  labs(
    y = "Bottom Temperature", 
    x = "Year", 
    color = "Season",
    fill = "Season", 
    linetype = "Trend",
    title = "")



temp1 

Seasonal bottom temperatures have increased in both spring and fall for the Northeast Shelf.

Looking at different areas within the region, Fall temperatures are increasing across the board, but seasonally spring temperatures in SNE and MAB are stable.

D. All Trends Figure

# Show plot
( 
    # Size Spectra
    (bmspectra_trend_p + 
       labs(title = "") + 
       theme(
         strip.text.y = element_blank(),
         plot.margin = margin(0,2,0,2),
     legend.position = "none")) |
     # Median Weight
     (wt_plot + 
        theme(
          strip.text.y = element_blank(),
          plot.margin = margin(0,2,0,2),
          legend.position = "none")) |
     # Bottom Temperature
     temp1 + 
      theme(
        plot.margin = margin(0,2,0,2)) 
    ) +
  plot_layout(guides = "collect") & theme(legend.title = element_text(hjust = 0.5),
                                          legend.title.position = "top")

E. Small/Large Fish Indices

NOTE: Percentile location of the size distribution are determined using DescTools::Quantile() which uses numlen_adj to weight everything by abundance, properly accounting for abundance.

The following figure shows changes in seasonal large and small fish indices within each region. The large and small size thresholds are set for each species independently using their length distributions as in the table above.

The axes on these are: \[LFI = \frac{\sum({N_{tjsk} > size~threshold_k})} {\sum{N_{tjs}}}\]

for each: \(year_t\), \(region_j\), \(season_s\) & \(species_k\).

# Get percentiles based on length for each species
# # Get 90th or 95th percentile as size threshold
# Do them separately
species_size_quants <- wigley_in %>% 
  mutate(area = "Northeast Shelf") %>% 
  unite("species_area", c(area, comname), sep = "X", ) %>% 
  split(.$species_area) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.05, 0.1, 0.5, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "species_area") %>% 
  separate(col = species_area, into = c("area", "comname"), sep = "X") %>% 
  mutate(area = factor(area, levels = area_levels_long_all)) %>% 
  arrange(area)


# # Show the table
# species_size_quants %>% 
#   mutate_if(is.numeric, round, 2) %>% 
#   head(n = 10) %>% 
#   gt() %>% 
#   tab_header(title = "Wigley Community - Body Mass Quantiles (kg)") %>% tab_spanner(label = "Percentile", columns = -area) %>% 
#   cols_label("area" = "")

# Assign large/small indicators
wigley_lfi <- wigley_in %>% 
  left_join(select(species_size_quants, comname, `10%`, `90%`)) %>%   
  mutate(
    small = ifelse(ind_weight_kg < `10%`, T, F),
    large = ifelse(ind_weight_kg > `90%`, T, F))


# Get total abundance in each area*year*season
total_abundbio <- bind_rows(
    mutate(wigley_lfi, area = "Northeast Shelf"),
    wigley_lfi) %>% 
  group_by(area, est_year, season) %>% 
  summarise(
    total_abund = sum(numlen_adj, na.rm = T),
    .groups = "drop")
  
  
# Get large/small fish  abundances
lf_abundbio <- bind_rows(
    mutate(wigley_lfi, area = "Northeast Shelf"),
    wigley_lfi) %>%
  group_by(area, comname, est_year, season) %>%
  summarise(
    small_abund = sum(numlen_adj * small),
    large_abund = sum(numlen_adj * large),
    .groups = "drop") %>% 
  group_by(area, est_year, season) %>% 
  summarise(
    small_abund = sum(small_abund, na.rm = T),
    large_abund = sum(large_abund, na.rm = T),
    .groups = "drop")

# Join and divide
lfi <- left_join(total_abundbio, lf_abundbio) %>%
  mutate(
    LFI = large_abund / total_abund,
    SFI = small_abund / total_abund)  %>%
  mutate(
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) 

# Plot them as loing-term average

# Plot the SFI
sfi_p <- lfi %>% 
  ggplot() +
  geom_point(
    aes(est_year, SFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, SFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(labels = label_percent()) +
  theme(strip.text.y = element_blank()) +
  facet_grid(area~., scales = "fixed", labeller = label_wrap_gen(width = 8)) +
  labs(
    x = "Year",
    y = "Percent of Total Abundance",
    title = "Wigley Species",
    subtitle = "Small Fish Index",
    fill = "Season",
    color = "Season",
    linetype = "")


# Plot the LFI
lfi_p <- ggplot(lfi) +
  geom_point(
    aes(est_year, LFI, color = season),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(est_year, LFI, color = season, linetype = "5-Year Moving Average"), 
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_fill_gmri() +
  scale_color_gmri() +
  scale_y_continuous(
    labels = label_percent(), 
    expand = expansion(add = c(0,0.05))) +
  facet_grid(area~., scales = "fixed", labeller = label_wrap_gen(width = 8)) +
  labs(
    x = "Year",
    y = "",
    subtitle = "Large Fish Index",
    fill = "Season",
    color = "Season",
    linetype = "")



# Plot together
(sfi_p | lfi_p) + plot_layout(guides = "collect")

From these figures we can see that there is a general pattern where larger individuals (relative to their conspecific size-distribution) make up a smaller fraction of seasonal abundances.

There is something strange going on with GB in the 2000’s where the LFI increases for a time, but this is no longer the case.

Upper Percentile Drift of the Size distribution

These metrics are an attempt to capture where movements at the edges of the size distribution: Is the size distribution becoming wider or more narrow?

There is a sentiment that the 90th percentile fish is smaller now than it has been in the past. This is true in some circumstances but not all. Indeed, the 90th percentile mass for fishes in the GOM was nearly 2x what it is today. A similar change in the upper percentile location is visible in Georges Bank as well. At the shelf-wide scale and in the southern regions this is less true.

10th percentile weights seem stable and are likely less informative due to gear selectivity preventing consistent sampling at the smallest sizes. I find it more probable that the size distribution changes are related to fluctuations in larger individuals than those at the lower end of the distribution.

# Do them separately
yearly_size_quants <- bind_rows(
    mutate(wigley_in, 
           area = "Northeast Shelf"),
    wigley_in) %>% 
  unite(col = "group_var", area, est_year, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.90, .95)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all)) %>% 
  arrange(area)


# Reshape to make filtering easier
quants_long <- yearly_size_quants %>% 
  pivot_longer(
    cols = -c(area, est_year), 
    names_to = "Percentile", 
    values_to = "biomass_kg") 


#-------------

# Take the data excluding dogfish
no_spiny <- filter(wigley_in, comname != "spiny dogfish")


# Do them separately
spineless_size_quants <- bind_rows(
    mutate(no_spiny, area = "Northeast Shelf"),
    no_spiny) %>% 
  unite(col = "group_var", area, est_year, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1,  0.5, 0.90)) %>% 
      t() %>% 
      as_tibble() 
    
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all)) %>% 
  arrange(area)


# Reshape to make filtering easier
spineless_quants_long <- spineless_size_quants %>% 
  pivot_longer(
    cols = -c(area, est_year), 
    names_to = "Percentile", 
    values_to = "biomass_kg")

# Separate axes
quants_10 <- quants_long %>% filter(Percentile == "10%") %>% mutate(Percentile = "10th Percentile")
spiny_10 <-  spineless_quants_long %>% filter(Percentile == "10%") %>% mutate(Percentile = "10th Percentile")
q10_p <- ggplot() +
  geom_point(
    data = quants_10,
     aes(est_year, biomass_kg*1000, 
         color = "Wigley Community"),
     size = 0.8, alpha = 0.5) +
  geom_ma(
    data = quants_10,
    aes(x = est_year, y = biomass_kg*1000, 
        color = "Wigley Community",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_point(
    data = spiny_10,
    aes(est_year, biomass_kg*1000, 
        color = "Spiny Dogfish Removed"),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = spiny_10,
    aes(x = est_year, y = biomass_kg*1000,
        color = "Spiny Dogfish Removed",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_color_grey() +
  scale_y_continuous(labels = label_number(suffix = "g")) +
  guides(
    linetype = guide_legend(title.position = "top", title.hjust = 0.5),
    color = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      nrow = 2)) +
  theme(
    strip.text.y = element_blank(),
    plot.margin = margin(0,0,0,0)) +
  facet_grid(area~Percentile, labeller = label_wrap_gen(width = 8)) +
  labs(
    y = "Bodymass", 
    x = "Year", 
    linetype = "Linetype",
    color = "Community")



# Separate axes for 90
quants_90 <- quants_long %>% filter(Percentile == "90%") %>% mutate(Percentile = "90th Percentile")
spiny_90 <-  spineless_quants_long %>% filter(Percentile == "90%") %>% mutate(Percentile = "90th Percentile")
q90_p <-  ggplot() +
  geom_point(
    data = quants_90,
     aes(est_year, biomass_kg, 
         color = "Wigley Community"),
     size = 0.8, alpha = 0.5) +
  geom_ma(
    data = quants_90,
    aes(x = est_year, y = biomass_kg, 
        color = "Wigley Community",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  geom_point(
    data = spiny_90,
    aes(est_year, biomass_kg, 
        color = "Spiny Dogfish Removed"),
    size = 0.8, alpha = 0.5) +
  geom_ma(
    data = spiny_90,
    aes(x = est_year, y = biomass_kg,
        color = "Spiny Dogfish Removed",
        linetype = "5-Year Moving Average"),
    n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  scale_color_grey() +
  scale_y_continuous(labels = label_number(suffix = "kg")) +
  guides(
    linetype = guide_legend(title.position = "top", title.hjust = 0.5),
    color = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      nrow = 2)) +
  theme(plot.margin = margin(0,0,0,0)) +
  facet_grid(area~Percentile, labeller = label_wrap_gen(width = 8)) +
  labs(
    y = "", 
    x = "Year", 
    linetype = "Linetype",
    color = "Community")



# Assemble
(q10_p | q90_p) + plot_layout(guides = "collect", axes = "collect_x") 

Dogfish & Seasonal Size Distributions

From exploratory analyses we know that spiny dogfish are a large fraction of the biomass in some seasonal communities. Because of this, their behavior has a large impact on the community size structure. However despite their outsized impact they are not the whole story and there is a danger in over-emphasizing their role as they are not strictly benthic and are known to migrate within and beyond the bounds of the survey sampling area.

By excluding spiny dogfish from the community the 2000’s era increase in 90th percentile size for the Gulf of Maine goes away. Possibly giving support to sentiments that spiny dogfish are occupying the body-size-niche perviously held by larger gadids.

# remake the above plot as seasonal

# Do them separately
wigley_season_qs <- bind_rows(
    mutate(wigley_in, 
           area = "Northeast Shelf"),
    wigley_in) %>% 
  unite(col = "group_var", area, est_year, season, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.5, 0.90)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year", "season"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) %>% 
  arrange(area) %>% 
  pivot_longer(
    cols = -c(area, est_year, season), 
    names_to = "Percentile", 
    values_to = "biomass_kg")  %>% 
  group_by(area, season, Percentile) %>% 
  arrange(est_year) %>% 
  mutate(kg_roll = zoo::rollapply(biomass_kg, 5, mean, na.rm = T, align = "right",  fill = NA))


#-------------

# Take the data excluding dogfish
no_spiny <- filter(wigley_in, comname != "spiny dogfish")


# Do them without dogfish separately
spiny_season_qs <- bind_rows(
    mutate(no_spiny, area = "Northeast Shelf"),
    no_spiny) %>% 
  unite(col = "group_var", area, est_year, season, sep = "X", remove = T) %>% 
  split(.$group_var) %>% 
  map_dfr(function(x,y){
    
    DescTools::Quantile(
      # x$length_cm, 
      x$ind_weight_kg,
      weights = x$numlen_adj, 
      probs = c(0.1, 0.5, 0.90)) %>% 
      t() %>% 
      as_tibble() 
  }, .id = "group_var") %>% 
  separate(group_var, into = c("area", "est_year", "season"), sep = "X") %>% 
  mutate(
    est_year = as.numeric(est_year),
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) %>% 
  arrange(area) %>% 
  pivot_longer(
    cols = -c(area, est_year, season), 
    names_to = "Percentile", 
    values_to = "biomass_kg") %>% 
  group_by(area, season, Percentile) %>% 
  arrange(est_year) %>% 
  mutate(kg_roll = zoo::rollapply(biomass_kg, 5, mean, na.rm = T, align = "right",  fill = NA))

Seasonal Size Percentiles

The seasonality in the 90th percentile shift with and w/o spiny dogfish does not speak to migration as much as I thought it might originally.

It is capturing a more general absence of dogfish (and other larger individuals) in MAB during the fall. The departure of larger individuals from the area seasonally is possibly related to the breakdown of the Mid-Atlantic Cold pool, which breaks down over the summer.

# Separate axes for 90
quants_90 <- wigley_season_qs %>% 
  filter(Percentile == "90%") %>% 
  mutate(Percentile = "90th Percentile")
spiny_90 <-  spiny_season_qs %>% 
  filter(Percentile == "90%") %>% 
  mutate(Percentile = "90th Percentile")


# 90th percentile dogfish changes
q90_df_p <-  ggplot() +
  geom_area(
    data = quants_90,
    aes(x = est_year, y = kg_roll, 
        fill = "Complete Community"),
    alpha = 0.6, color = "gray70") +
  geom_area(
    data = spiny_90,
    aes(x = est_year, y = kg_roll, 
        fill = "Spiny Dogfish Removed"),
    alpha = 0.6, color = "gray20") +
  geom_image(
    data = filter(
      quants_90, 
      est_year == 2008, 
      !(area == "Mid-Atlantic Bight" & season == "Fall")), 
    aes(est_year, kg_roll*1.65, 
        image = here::here("Data/phylopics/squalus_phylopic.png")),
    size = 0.35) +
  scale_fill_manual(values = c("gray70", "gray20")) +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  scale_y_continuous(labels = label_number(suffix = "kg")) +
  #coord_flip() +
  guides(
    linetype = guide_legend(title.position = "top", title.hjust = 0.5),
    color = guide_legend(
      title.position = "top",
      title.hjust = 0.5,
      nrow = 2)) +
  theme(plot.margin = margin(0,0,0,0),
        strip.background = element_rect(color = "white", linewidth = 0.2)) +
  facet_nested(
    area~Percentile*season,
    labeller = labeller(
      area = label_wrap_gen(width = 8),
      Percentile = ~str_c(.x, " of Body Mass Distribution"))) +
  labs(
    y = "Individual Bodymass", 
    x = "Year", 
    fill = "Community")
q90_df_p

The above figure shows that the community distribution without dogfish widens in the spring and narrows in the fall.

Percent Dogfish

It is maybe worth noting just how much of the community biomass is held in spiny dogfish, which is frequently above 50%.

dogfish_fraction <- bind_rows(
    mutate(
      wigley_in, 
      area = "Northeast Shelf"),
    wigley_in) %>% 
  group_by(area, est_year, season) %>% 
  summarise(
    dogfish_abund = sum(numlen_adj * (comname == "spiny dogfish"), na.rm = T),
    total_abund   = sum(numlen_adj, na.rm = T),
    dogfish_bio   = sum(sum_weight_kg * (comname == "spiny dogfish"), na.rm = T),
    total_bio     = sum(sum_weight_kg, na.rm = T),
    .groups = "drop") %>% 
  mutate(
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")),
    dogfish_abund_fraction = (dogfish_abund / total_abund),
    dogfish_bio_fraction = (dogfish_bio / total_bio))

dogfish_fraction %>% 
  pivot_longer(cols = ends_with("fraction"), names_to = "metric", values_to = "fraction") %>% 
  mutate(metric = if_else(str_detect(metric, "abund"), "Abundance", "Biomass")) %>% 
ggplot(aes(est_year, fraction)) +
  geom_area(aes(color = season, fill = season), alpha = 0.2) +
  #geom_line(aes(color = season, fill = season), alpha = 0.2) +
  facet_nested(area~metric*season) +
  scale_y_continuous(
    labels = label_percent(),
    limits = c(0,1), expand = expansion(add = c(0,0))) +
  scale_color_gmri() +
  scale_fill_gmri() +
  labs(y = "Community Fraction",
       x = "Year",
       title= " Dogfish as a Fraction of Total Abundance/Bio",
       color = "Season",
       fill = "Season")

Bottom Temperature CCF

Shelf-Wide:

Spring bottom temperatures lead size distributions by 4 or 3 years. The spring temperature part may be important because that is the annual low i.e. the baseline is shifting.

Regionally:

Bottom temperature appears fairly centered (sequential significant correlations at or around 0 lag/lead) in GOM & GB where it is more highly correlated and there is more consistent correlation direction and magnitude at sequential lags.

For MAB-Spring a 1 year of bottom temperature may be important.

In the Fall the GOM has its largest negative correlation centered at a 2-year lead.

# Function to grab the correlation data and lag data
get_ccf_vector <- function(x,y, lagmax = 10){
  
  # Run the ccf
  ccf_data <- ccf(x, y, plot= F , lag.max = lagmax)
  
  # Get the signif:
  # https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/
  # Not 100% sure n is the same for ccf as it is for acf, but...
  significance_level <- qnorm((1 + 0.95)/2)/sqrt(sum(!is.na(x)))
  
  data.frame(
    "acf"    = ccf_data$acf,
    "lag"    = ccf_data$lag,
    "sigpos" = significance_level,
    "signeg" = significance_level*-1
  )
}


# Run Local CCF on Bottom Temperature
btemp_ccf <- wigley_bmspectra_model_df %>% 
  unite(col = "regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(~get_ccf_vector(
    x = .x$bot_temp, 
    y = .x$b), .id = "regseas") %>% 
  separate(col = "regseas", into = c("area", "season"), sep = "XX") %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))



ggplot(btemp_ccf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(btemp_ccf, sig_flag), color = "black") +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Bottom Temperature CCF",
    x = "Bottom Temperature Lag (years)", 
    y = "ACF", )

Landings CCF

Its plausible that landings have a lagged effect so here are the CCFs for that.

These visually show a pattern of sustained autocorrelation, but in most instances the correlation to landings is such that the slope of the size distribution changes leads the landings AND that higher landings are positively correlated with shallow spectra.

# Run Local CCF on Bottom Temperature
landings_ccf <- wigley_bmspectra_model_df %>% 
  drop_na(total_weight_lb) %>% 
  unite(col = "regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(~get_ccf_vector(
    x = .x$total_weight_lb, 
    y = .x$b), .id = "regseas") %>% 
  separate(col = "regseas", into = c("area", "season"), sep = "XX") %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall")))



ggplot(landings_ccf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(landings_ccf, sig_flag), color = "black") +
  scale_x_continuous(expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Total Landings CCF",
    x = "Total Landings Lag (years)", 
    y = "ACF", )

Spectra Slope ACF

Just as an important check, what is the temporal autocorrelation structure in spectra slopes?

There is a significant positive correlation at 1 year lead in 4/10 regional timeseries. I think the nature of the data is enough to justify adding the AR1() error structure, and this is confirms its inclusion over a year random effect.

There is also six-year lead in 3 of them which is curious, and could be explored as to what is going on there with a 6-year periodicity.

Then lastly the Fall Gulf of Maine spectra has significant correlations up to three years out.

I am of two minds about this:

1. This could be a clue to some state change, as increased autocorrelation is sometimes used in the literature as a signal of hysteresis.

2. OR This could be evidence that Gulf of Maine’s fall community is the exception, and is the ONLY community sampled sufficiently/consistently.

We are in a sense (hopefully) conducting a repeated measure of some community in time, and so there should hopefully be some sort of autocorrelation. The absence of autocorrelation in these communities could

# Pull out acf stuff
spectra_acf <- wigley_bmspectra_model_df %>% 
  unite("regseas", area, season, sep = "XX") %>% 
  split(.$regseas) %>% 
  map_dfr(function(x){
    # drop NA
    x <- arrange(x, est_year) %>% drop_na(b)
    
    # Do ACF
    x_acf <- acf(x$b, plot = F, lag.max = 10)
    
    # Get the signif:
    # https://www.squaregoldfish.co.uk/programming/r_acf_significance.md/
    # Not 100% sure n is the same for ccf as it is for acf, but...
    significance_level <- qnorm((1 + 0.95)/2)/sqrt(sum(!is.na(x$b)))
    
    data.frame(
      "acf"    = x_acf$acf,
      "lag"    = x_acf$lag,
      "sigpos" = significance_level,
      "signeg" = significance_level*-1
    )
    
  }, .id = "regseas") %>% 
  separate("regseas", into = c("area", "season"), sep = "XX")  %>% 
  mutate(
    # Flag when it crosses threshold
    sig_flag = ifelse(acf < signeg | acf > sigpos, T, F),
    # Set Factor Levels
    area = factor(area, levels = area_levels_long_all),
    season = factor(season, levels = c("Spring", "Fall"))) 



# Plot things
ggplot(spectra_acf, aes(lag, acf)) +
  geom_hline(yintercept = 0, color = "gray25", linewidth = 1) +
  geom_vline(xintercept = 0, color = "gray25", linewidth = 1) +
  geom_ribbon(
    aes(ymin = signeg, ymax = sigpos), 
        alpha = 0.2, linetype = 2, fill = "lightgray",
    color = "gray65") +
  geom_col(alpha = 0.65) +
  geom_col(
    data = filter(spectra_acf, sig_flag), color = "black") +
  scale_x_continuous(
    limits = c(-1,10),
    breaks = c(1:9),
    expand = expansion(add = c(0,0))) +
  facet_grid(area~season) +
  labs(
    title = "Spectra Slope ACF",
    x = "Lag (years)", 
    y = "ACF", )

Spring and Fall Correlated ??

I can’t believe I never checked this, but are spring and fall communities correlated within an area? My sense is that they should be, and it might be meaningful to point out if they are not.

Only one region (GOM) has “memory” of the community six months before, and has a significant (but small) correlation between spring and fall communities.

Taken with the ACF info, there is consistency that the Gulf of Maine spectra exhibits more sign of temporal autocorrelation and the spectra between season are correlated.

This does however add doubt to how useful these other indices are, in the sense that they may not be capturing change of a community over time.

wigley_bmspectra_model_df %>% 
  select(est_year, area, season, b) %>% 
  pivot_wider(names_from = "season", values_from = "b") %>% 
  ggplot(aes(Spring, Fall)) +
  geom_point() +
  geom_smooth(method = "lm") +
  ggpmisc::stat_poly_eq(ggpmisc::use_label(c("P", "R2")), formula = y ~ x) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.5))) +
  facet_wrap(~area, ncol = 2, scales = "free") +
  labs(title = "Spring-Fall Spectra Correlation")

Summary of Trends

Size spectra steepening and median size declines are fairly consistent with one another. Long-term trends are only present at regional scales, and not ocurring shelf-wide for this community of 74 species.

There is something unique happening during the spring in the Mid-Atlantic Bight, where median size is growing and the size spectra is shallowing. This region exhibits large seasonal differences in community composition and large fluctuations in bottom temperature partly due to the mid-atlantic cold-pool.

From the large fish indices I see evidence that individuals within each species are smaller now than they had been, and that larger individuals are rarer in the overall community.

The shift in the 90th percentile size also supports this idea of smaller individuals being more common. The differences seen when excluding spiny dogfish gives some additional evidence that this is particularly true for non-dogfish species. This may also support the idea that dogfish species replacement has ocurred and spiny dogfish are now occupying that larger size class niche space.

From the temperature CCf plot my takeaways are that spring bottom temperature matters more for the shelf-wide size distribution and may operate with a lag of several years. Regionally, the lag seems less relevant and a negative correlation is more evident in GOM+GB.

From the landings CCF I get the sense that total landings aren’t a major driver of slope steepening in our area, or they operate at very long lead times 8-10 years. Many of the significant correlations are in the opposite direction of how we think landings would impact the size distribution (removal of larger individuals leads to steepening), AND they are such that they follow the size distribution slope rather than lead it.

---

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

Temperature and Size Structure

There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.

• Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
• James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes

The metabolic theory of ecology provides a basis for “why” the above rules might be observed.

• Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.

Seasonal Movements and Metabolic Efficiency: MTE details how metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.

This theory would suggest that larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and to reduce the energetic costs of thermoregulation in warmer periods.

From what I see at first-glance the size distribution along the shelf is responding consistently with temperature changes (Spring to Fall) on shorter time scales than hypotheses that predict temperature related changes in growth.

wigley_bmspectra_model_df %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b, color = season)) +
  geom_point(aes(shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  ggforce::geom_mark_ellipse(aes(fill = season), alpha = 0.1) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  guides(
    shape = guide_legend(nrow = 2),
    fill = guide_legend(nrow = 2),
    color = guide_legend(nrow = 2)) +
  theme(legend.box = "horizontal") +
  labs(
    title = "Bergmann's Rule + MTE",
    subtitle = "The broad relationship between temperature and size distribution is bolstered by\nseasonal community changes in alignment with seasonal temperature differences.",
    shape = "Survey Area",
    fill = "Season",
    color = "Season",
    x = "Bottom Temperature")

Importance of Seasonality

The extent that the regional spectra trends follow a negative linear relationship to increasing temperature is largely due to seasonal differences.

We know that individuals may migrate and/or occupy different habitats seasonally, leading to situations where the community compositions from spring to fall are an apples to oranges comparison. Seasonal movements may be be partially explained by expectations from the MTE, but they also could be separate ecological processes like seasonal foraging or spawning opportunities.

Either way, seasonal differences are large and they justify treating survey seasons as separate groups with independent trends.

Anomalies Model

I wanted to try and avoid conflating that a relationship between spectra and bottom temperature means that “warming” has altered the size spectra. The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.

This is the model I would suggest gets the most at “Warming” impacting regional size distributions:

Fixed effects:

• Season
• Survey area
• Bottom temperature anomaly

Error structure: AR1 autocorrelative errors

By using the season * region interaction groups, we can characterize the typical bottom temperature for each combinations using the long-term seasonal average, and quantify warming using departures from those long-term averages.

The model looks like this as an equation (90% sure), where \(\alpha\) is the intercept, \(\beta x_{tjs}\) is our coefficient for a continuous independent variable (ex. bottom temperature), and \(e_{tjs\)` is our error:

\[y_{tjs} = \alpha + \beta x_{tjs} + e_{tjs}\] and our error depends on the previous timestep:

\[e_{tjs} = \phi e_{t-1,js} + w_{tjs}\] With \(w_{tjs}\) adding noise.

\[w_{tjs} \sim N(0, \sigma^2) \] All of this happening for each \(year_t\) within \(region_j\) & \(season_s\).

I need to consult the Zuur paper about equations in ecology again, but this seems close.

# Just drop the random effect part for now
anom_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp_anom), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(anom_model_ar)$corStruct


# # check the model
# check_model(anom_model_ar)
# plot(check_collinearity(anom_model_ar))
# plot(
#   check_collinearity(
#     nlme::gls(
#       b ~ area + season + scale(bot_temp_anom), #+ scale(log(total_weight_lb)), 
#       data = wigley_bmspectra_model_df, 
#       correlation = corAR1(form = ~ est_year | area/season)
# ))
#   )


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    anom_model_ar, 
    terms = ~ bot_temp_anom * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = anom_model_ar,
  specs =  ~ area * season,
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop") %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
local_btemp_anoms <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Size Distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Based on this model I would say there is evidence of a warming effect in:
1. GOM-Spring
2. GB-Spring
3. GB-Fall
4. SNE-Fall

Bottom Temperature Model

We were originally using absolute temperatures, which are scaled, and I was curious with the AR error structure how different this would be than the anomaly model.

I just had to see if the same model is better performing or leads to different results when it knows absolute temperatures.

# Just drop the random effect part for now
temp_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct
# 
# # check the model
# check_model(temp_model_ar)

# # Collinearity without interactions
# plot(
#   check_collinearity(
#      nlme::gls(
#         b ~ area + season + scale(bot_temp), 
#         data = wigley_bmspectra_model_df, 
#         correlation = corAR1(form = ~ est_year | area/season))
# ))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    temp_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = temp_model_ar,
  specs =  ~ area * season,
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")  %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
local_btemp_results <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Size Distribution Exponent",
       x = "Seasonal Bottom Temperature")
local_btemp_results

Verdict: Temperature itself performs better, and is simpler to produce/explain.

data.frame(
  "Model" = c("Temperature Model", "Anomalies Model"),
  "AIC" = c(AIC(temp_model_ar), AIC(anom_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()

Model	AIC	Delta-AIC
Temperature Model	-283.6706	0.00000
Anomalies Model	-261.7763	21.89431

Thoughts:

Consistent seasonal shifts in the Shelf-wide spectra suggests to me that the community is more responsive to within-year time scale factors and maybe less sensistive to long-term trends in temperature than I originally expected. This also raises questions about how the spring/fall size distribution differences are being achieved and whether to view them as expressions of the same “community” or separately.

From the decadal variability work we have reason to believe that individual species are both not (perfectly) tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

In light of that knowledge, there must be some compensatory mechanisms happening across this broader community for the temperature+spectra relationship to be so limited/situational.

Region * Season Pairwise Posthoc

The following figures and tables compare marginal means and post-hoc comparisons for the model above using bottom temperature.

What are the differences across spatial scales and seasonally?

# Region and seasonal posthocs
regseas_phoc <- emmeans(
  temp_model_ar, 
  list(pairwise ~ area + season), 
  adjust = "tukey",
  type = "response")



# Custom Plot
(intercepts_emmeans <- regseas_phoc$`emmeans of area, season` %>% 
    as_tibble() %>%
    ggplot() +
    geom_pointrange(
      aes(area, emmean, ymin = lower.CL, ymax = upper.CL, color = season),
      position = position_dodge(width = 0.25), alpha = 1) +
    #coord_flip() +
    scale_color_gmri() +
    labs(
      y = "Marginal Mean: Mass Spectra Slope",
      x = "Area",
      title = "Wigley Species",
      color = "Season",
      subtitle = "Group marginal means plot"))

# Table for pairwise test results
regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  arrange(p.value) %>% 
  gt() %>% 
  tab_header(title = "Factor Groups Pairwise Post-Hoc") %>% 
  tab_style(
    style = list(
      cell_fill(color = "#9AFF9A"), 
      cell_text(color = "black")),
    locations = cells_body(
      rows = p.value <= 0.05) 
  )

Factor Groups Pairwise Post-Hoc
Contrast	estimate	SE	df	t.ratio	p.value
Northeast Shelf Spring - (Mid-Atlantic Bight Fall)	0.943	0.130	471.630	7.252	0.000
Gulf of Maine Spring - (Mid-Atlantic Bight Fall)	0.878	0.134	475.331	6.549	0.000
Georges Bank Spring - (Mid-Atlantic Bight Fall)	0.862	0.150	474.065	5.731	0.000
Southern New England Spring - (Mid-Atlantic Bight Fall)	0.931	0.119	474.694	7.805	0.000
(Mid-Atlantic Bight Spring) - (Mid-Atlantic Bight Fall)	0.978	0.109	473.235	9.014	0.000
Northeast Shelf Fall - Georges Bank Fall	-0.384	0.082	474.790	-4.660	0.000
Northeast Shelf Fall - (Mid-Atlantic Bight Fall)	0.733	0.117	439.871	6.287	0.000
Gulf of Maine Fall - (Mid-Atlantic Bight Fall)	1.025	0.118	442.826	8.695	0.000
Georges Bank Fall - (Mid-Atlantic Bight Fall)	1.117	0.122	477.971	9.173	0.000
Southern New England Fall - (Mid-Atlantic Bight Fall)	0.742	0.125	476.724	5.913	0.000
(Mid-Atlantic Bight Spring) - Northeast Shelf Fall	0.245	0.061	458.918	4.011	0.003
Georges Bank Fall - Southern New England Fall	0.375	0.094	475.495	3.976	0.003
Northeast Shelf Fall - Gulf of Maine Fall	-0.292	0.077	460.546	-3.815	0.006
Gulf of Maine Fall - Southern New England Fall	0.283	0.089	471.001	3.173	0.051
(Mid-Atlantic Bight Spring) - Southern New England Fall	0.237	0.077	471.559	3.091	0.065
Southern New England Spring - Northeast Shelf Fall	0.198	0.079	470.459	2.518	0.262
Gulf of Maine Spring - Georges Bank Fall	-0.239	0.106	477.157	-2.261	0.417
Northeast Shelf Spring - Northeast Shelf Fall	0.209	0.094	463.762	2.226	0.441
Southern New England Spring - Georges Bank Fall	-0.186	0.086	477.144	-2.156	0.489
Southern New England Spring - Southern New England Fall	0.190	0.091	474.120	2.079	0.543
Georges Bank Spring - Georges Bank Fall	-0.255	0.126	476.856	-2.031	0.578
(Mid-Atlantic Bight Spring) - Georges Bank Fall	-0.139	0.070	476.051	-1.970	0.621
Northeast Shelf Spring - Southern New England Fall	0.201	0.105	469.083	1.917	0.657
Northeast Shelf Spring - Georges Bank Fall	-0.174	0.100	471.737	-1.739	0.773
Gulf of Maine Spring - Northeast Shelf Fall	0.145	0.100	473.576	1.454	0.909
Gulf of Maine Spring - Gulf of Maine Fall	-0.147	0.101	476.177	-1.453	0.909
Georges Bank Spring - Gulf of Maine Fall	-0.163	0.122	476.375	-1.340	0.944
Gulf of Maine Spring - Southern New England Fall	0.136	0.110	471.249	1.243	0.965
Southern New England Spring - Gulf of Maine Fall	-0.094	0.081	477.748	-1.165	0.977
Gulf of Maine Spring - (Mid-Atlantic Bight Spring)	-0.100	0.090	473.850	-1.112	0.983
Gulf of Maine Fall - Georges Bank Fall	-0.092	0.084	444.940	-1.092	0.985
Georges Bank Spring - Northeast Shelf Fall	0.129	0.121	469.365	1.066	0.988
Georges Bank Spring - (Mid-Atlantic Bight Spring)	-0.116	0.113	477.514	-1.032	0.990
Georges Bank Spring - Southern New England Fall	0.120	0.129	431.576	0.930	0.995
Northeast Shelf Spring - Gulf of Maine Fall	-0.083	0.096	469.882	-0.864	0.997
Southern New England Spring - (Mid-Atlantic Bight Spring)	-0.047	0.066	477.952	-0.711	0.999
(Mid-Atlantic Bight Spring) - Gulf of Maine Fall	-0.047	0.063	470.330	-0.738	0.999
Northeast Shelf Spring - Gulf of Maine Spring	0.064	0.115	477.403	0.560	1.000
Northeast Shelf Spring - Georges Bank Spring	0.081	0.134	475.084	0.604	1.000
Northeast Shelf Spring - Southern New England Spring	0.011	0.097	408.990	0.115	1.000
Northeast Shelf Spring - (Mid-Atlantic Bight Spring)	-0.036	0.084	474.661	-0.427	1.000
Gulf of Maine Spring - Georges Bank Spring	0.016	0.138	433.350	0.119	1.000
Gulf of Maine Spring - Southern New England Spring	-0.053	0.103	477.368	-0.517	1.000
Georges Bank Spring - Southern New England Spring	-0.070	0.123	476.681	-0.564	1.000
Northeast Shelf Fall - Southern New England Fall	-0.008	0.088	417.718	-0.096	1.000

regseas_phoc %>% 
  pluck(2) %>% 
  as.data.frame() %>% 
  rename(Contrast = `1`) %>% 
  mutate_if(is.numeric, ~round(.x, 3)) %>% 
  filter(p.value < 0.05)  %>% 
  separate(Contrast, into = c("group1", "group2"), sep = " - ") %>% 
  mutate(
    group1 = str_remove_all(group1, "[(]|[)]"),
    group2 = str_remove_all(group2, "[(]|[)]"),
    statement = 
      if_else(estimate > 0, " is shallower than ", " is steeper than ")) %>% 
  select(group1, statement, group2) %>% 
  arrange(statement, group1) %>% 
  gt()

group1	statement	group2
Georges Bank Fall	is shallower than	Southern New England Fall
Georges Bank Fall	is shallower than	Mid-Atlantic Bight Fall
Georges Bank Spring	is shallower than	Mid-Atlantic Bight Fall
Gulf of Maine Fall	is shallower than	Mid-Atlantic Bight Fall
Gulf of Maine Spring	is shallower than	Mid-Atlantic Bight Fall
Mid-Atlantic Bight Spring	is shallower than	Northeast Shelf Fall
Mid-Atlantic Bight Spring	is shallower than	Mid-Atlantic Bight Fall
Northeast Shelf Fall	is shallower than	Mid-Atlantic Bight Fall
Northeast Shelf Spring	is shallower than	Mid-Atlantic Bight Fall
Southern New England Fall	is shallower than	Mid-Atlantic Bight Fall
Southern New England Spring	is shallower than	Mid-Atlantic Bight Fall
Northeast Shelf Fall	is steeper than	Gulf of Maine Fall
Northeast Shelf Fall	is steeper than	Georges Bank Fall

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

# Plot annual totals
wigley_bmspectra_model_df %>% 
  filter(season == "Spring") %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.8) +
  geom_ma(
    aes(color = "5-Year Moving Average"), 
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries", linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Landings Autoregressive model

I wanted to follow on the bottom temperature anomalies model here for landings:

# Just drop the random effect part for now
f_model_ar <- nlme::gls(
  b ~ area * season + area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
    dplyr::filter(est_year > 1978),
  correlation = corAR1(form = ~ est_year | area/season)
)



# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    f_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = f_model_ar,
  ~area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite")


# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = f_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")  %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
landings_ar_plot <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(mod2_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free") +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

# temp model is better
# AIC(temp_model_ar)
# AIC(f_model_ar)

This is our old friend “higher landings when spectra was shallower” that we’ve seen before.

Temperature and Landings Together

If we want to just throw them both in because we’ve done what we can and want to hold them both in the frame at once we get this:

# Just drop the random effect part for now
full_model_ar <- nlme::gls(
  # Model
  model = b ~ area * season * scale(bot_temp) + area * scale(log(total_weight_lb)), 
  
  # Data
  data = wigley_bmspectra_model_df %>%
    drop_na(total_weight_lb) %>% 
    dplyr::filter(area != "Northeast Shelf"),
  
  # Autocorrelation
  correlation = corAR1(form = ~ est_year | area/season))

# # Collinearity without interactions
# plot(
#   check_collinearity(
#     nlme::gls(
#         b ~ area + season + scale(bot_temp) +  scale(log_land),
#         data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>%
#                  filter(area != "Northeast Shelf") %>% 
#           mutate(log_land = log(total_weight_lb)),
#         correlation = corAR1(form = ~ est_year | survey_area/season)
#       )
#     )
#   )


# # confidence intervals on phi
# intervals(full_model_ar)$corStruct


# temp model is much better
# AIC(temp_model_ar)
# AIC(full_model_ar)

Before proceeding it is worth noting that this model is less parsimonious than a temperature or anomaly only model:

data.frame(
  "Model" = c(
    "Temperature Model", 
    "Anomalies Model", 
    "Landings Model", 
    "Landings and Temperature"), 
  "AIC" = c(
    AIC(temp_model_ar), 
    AIC(anom_model_ar),
    AIC(f_model_ar),
    AIC(full_model_ar))) %>% 
  mutate(`Delta-AIC` = AIC - min(AIC)) %>% 
  gt()

Model	AIC	Delta-AIC
Temperature Model	-283.6706	0.00000
Anomalies Model	-261.7763	21.89431
Landings Model	-230.3029	53.36772
Landings and Temperature	-159.7297	123.94097

# Clean up the plot:
# Plot the predictions over data
temp_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
temp_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area * season,
  mode = "appx-satterthwaite",
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  filter(area != "Northeast Shelf") %>% 
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")  %>% 
  mutate(season = factor(season, levels = c("Spring", "Fall")))



# Plot over observed data
local_btemp <- temp_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(temp_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area!= "Northeast Shelf"),
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Size Distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp

# Clean up the plot:
# Plot the predictions over data
f_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
f_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite"
  ) %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- f_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(f_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area != "Northeast Shelf"),
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()

Cold Pool Index

There are several cold pool indices available that might be helpful in adding context to changes in MAB. Below are the Dupontavice et al. 2022 methods that are used in the SOE report.

The Cold Pool Index (Model_CPI) was adapted from Miller et al. (2016). Residual temperature was calculated in each grid cell, \(i\), in the Cold Pool domain as the difference between the average bottom temperature at the year \(y\) (\(T_{i,y}\)) and the average bottom temperature over the period 1972–2019 (\(\bar{T}_{i, 1972-2019}\)⁠⁠) between June and September. \(CPI\) was calculated as the mean residual temperature over the Cold Pool domain such that

\[CPI_y~=~\frac{\sum_{i=1}^{n}~(T_{i,y}-\bar{T}_{i,~1972-2019})}{n}\]

where n is the number of grid cells over the Cold Pool domain.

Takeaway:

\(CPI\) was calculated as the mean residual temperature over the Cold Pool domain

From this index we can see that the temperatures over the cold pool domain (which include the MAB) have increased over time.

cpool <- ecodata::cold_pool
cpool %>% 
  filter(!str_detect(Var, "se_")) %>% 
  ggplot(aes(Time, Value)) +
  geom_line() +
  facet_grid(Var~EPU, scales = "free")

# # There's an SF
# cpool_sf <- ecodata::cold_pool_sf %>% 
#   mutate(id = row_number())
# 
# ggplot(cpool_sf) + geom_sf(aes(fill = id), alpha = 0.2)

---

Driver Correlations